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ABSTRACT 

We study the dynamics of a planet on an orbit inclined with respect to a 
disc. If the initial inclination of the orbit is larger than some critical value, 
the gravitational force exerted by the disc on the planet leads to a Kozai 
cycle in which the eccentricity of the orbit is pumped up to large values and 
oscillates with time in antiphase with the inclination. On the other hand, 
both the inclination and the eccentricity are damped by the frictional force 
that the planet is subject to when it crosses the disc. We show that, by 
maintaining either the inclination or the eccentricity at large values, the Kozai 
effect provides a way of delaying alignment with the disc and circularization 
of the orbit. We find the critical value to be characteristically as small as 
about 20 degrees. Typically, Neptune or lower mass planets would remain on 
inclined and eccentric orbits over the disc lifetime, whereas orbits of Jupiter or 
higher mass planets would align and circularize. This could play a significant 
role in planet formation scenarios. 
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2 J. Teyssandier, C. Terquem and J. Papaloizou 
1 INTRODUCTION 

At the time of writing, more than 800 extrasolar planets have been detected around main 
sequence stars. Most of these objects have been observed through radial velocity measure- 
ments, and about one third of them have also been detected through transit measurements. 
In addition, the Kepler mission has so far released about 2300 planet candidates (Borucki 
et al. 2011, Batalha et al. 2012). Only a few of these objects have been confirmed by radial 
velocity measurements so far, but the number of false positive is expected to be very small. 
The rate at which extrasolar planets are being discovered is accordingly increasing very 
sharply. 

The angle between the sky projections of the stellar spin axis and the orbit normal 
(that we will call hereafter the projected inclination angle) has been measured using the 
Rossiter-McLaughlin effect for 53 planets (Albrecht et al. 2012 and references therein). 
About one third display significant misalignment with retrograde orbits being indicated in 
some cases. Misalignment, therefore, is common, at least for short period systems for which 
measurements have been made. Note however that the planets in this sample are rather 
massive. The lightest one has a mass of about 25 earth masses, whereas the other planets 
in the sample have masses ranging from a few tenths of a Jupiter mass to several Jupiter 
masses. By performing duration ratio statistics on the Kepler planetary candidates, Fabrycky 
et al. (2012) have concluded that pairs of planets in their sample are aligned to within a few 
degrees. As the Kepler catalog contains a majority of super-earth and Neptune-like planets, 
it is possible that misalignment is less common for lower mass planets. 

According to the commonly accepted core accretion scenario, in which planets form in 
a disc through the accretion of solid material into a core followed, in the case of giant 
planets, by the capture of a gaseous atmosphere, the orbit of planets should lie in the disc, 
and therefore in the equatorial plane of the star. In this model, migration due to tidal 
interaction between disc and planets can be invoked to explain the presence of giant planets 
on very short orbits. However, such an interaction would not account for the occurrence of 
any inclination of the orbit. Therefore, this scenario is not likely to apply to the hot Jupiter s 
that are observed to be on inclined orbits. Planet formation through fragmentation of the 
disc might be considered, but in that case the orbits are also expected to lie in the plane of 
the disc. 

One scenario that has been proposed for misaligning the orbits of hot Jupiters that have 
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formed in a disc lying in the stellar equatorial plane relies on gravitational interaction with 
a distant stellar or planetary companion that takes place after the disc has dissipated. When 
the two orbits are not coplanar, the secular perturbation exerted by the companion produces 
a Kozai cycle in which the eccentricity and the inclination of the inner planet's orbit vary in 
antiphase. If the pericenter distance gets small enough, the orbit may become circularized 
because of the tidal interaction with the central star, and the orbit shrinks while, in some 
circumstances, keeping a high inclination (Fabrycky & Tremaine 2007, Wu et al. 2007, Naoz 
et al. 2011). In this model, a distant companion is needed which, however, may not always 
be present in reality. Other models rely on planet-planet scattering (Chatterjee et al. 2008) 
or secular chaos (Wu & Lithwick 2011). 

It has been proposed that the disc itself may be misaligned with respect to the stellar 
equatorial plane. That could happen if the disc, at later times, were accreting material with 
angular momentum misaligned with that of the star, as considered by Bate et al. (2010). If at 
that stage there was enough mass in the disc to form planets, their orbits would naturally be 
inclined with respect to the stellar equatorial plane. A similar scenario was studied by Thies 
et al. (2011), who pointed out that close encounters of a disc accreting from an extended 
envelope with another star could result in the disc plane becoming tilted, possibly even 
to a retrograde orientation with respect to its original one. Planetary orbits inclined with 
respect to the stellar equatorial plane would also result if the stellar spin axis were tilted 
due to interaction with the disc (Lai et al. 2011, Foucart & Lai 2011). Note however that, 
by comparing stellar rotation axis inclination angles with the geometrically measured disc 
inclinations for a sample of eight debri-discs, Watson et al. (2011) have seen no evidence for 
a misalignment between the two. 

Another possibility would be that the misaligned planets have formed out of the disc 
through a fragmentation process occurring in the protostellar envelope while it collapses onto 
the protostar, as envisioned by Papaloizou & Terquem (2001) and Terquem & Papaloizou 
(2002). In this scenario, a population of planetary objects form rapidly enough that their 
orbits can undergo dynamical relaxation on a timescale on the order of a few 10 4 years. 
During the relaxation process, most of the objects are ejected, while one to a few planets 
become more bound to the central star. Formation of hot Jupiters through tidal circulariza- 
tion by the central star of a highly eccentric orbit may occur (Papaloizou & Terquem 2001). 
The orbits of the planets that are left in the system at the end of the relaxation process 
display a range of eccentricities and inclinations (Adams & Laughlin 2003). In this context, 
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the misaligned planets interact with the disc that has formed around the star during the 
protostellar envelope collapse. It is the subsequent dynamics of a system of this type, con- 
sisting of a gaseous disc together with a planet with an orbital plane misaligned with its 
midplane that we propose to investigate here. Note that we do not expect low mass planets 
to form according to the scenario of Papaloizou & Terquem (2001). On the other hand, 
in the scenario proposed by Bate et al. (2010), it may happen that subsequent discs with 
different inclinations form and dissipate, so that planets formed early in a disc may interact 
at later time with a disc with different orientation. 

In the present paper, we focus on a system with only one planet. Multiple systems will be 
studied in future publications. Some preliminary work on the interaction of a planet on an 
inclined orbit with a disc has been carried out by Terquem & Ajmia (2010). It was found that 
the gravitational force exerted by the disc onto the planet leads to a Kozai cycle in which 
the eccentricity and the inclination of the orbit vary periodically with large amplitudes. This 
indicates that a planet's orbit which is inclined to start with may achieve high eccentricity. In 
their calculations, Terquem & Ajmia (2010) adopted a two dimensional flat disc model and 
ignored the frictional force felt by the planet as it passes through the disc. In the present 
paper, we extend this work by modelling the disc in three dimensions and including the 
frictional force. Our goal is to understand under what circumstances inclined orbits can be 
maintained when the disc is present. We consider planets with masses ranging from Neptune 
mass to several Jupiter masses. The plan of the paper is as follows: 

In section I2.1[ we describe the three dimensional disc model used in the numerical sim- 
ulations and the calculation of the gravitational force exerted by the disc on a planet in 
an inclined orbit. A computationally convenient formulation, in terms of elliptic integrals, 
is presented in an appendix. We go on to give a brief review of the Kozai mechanism in 
section 12.21 In section 12.31 we give an expression for the frictional force exerted by the disc 
on the planet, and derive a damping timescale based on a simplified analysis in section EU 
In section [3j we present the results of numerical simulations of the interaction between a 
planet on an inclined orbit and a disc. We first perform simulations without the frictional 
force in section I3TT1 and compare the results to those of Terquem & Ajmia (2010) that were 
obtained for a two dimensional flat disc model. In section 13.21 we include friction. We show 
that, when the planet's orbit starts with a low inclination (less than about 23°) with re- 
spect to the disc, it becomes aligned with the disc plane and circularized by the frictional 
force. However, when the inclination is initially high enough, the Kozai effect is present. 
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This pumps up the eccentricity of the planet's orbit maintaining either the inclination or 
the eccentricity at large values. As a consequence, alignment of the orbital and disc planes, 
as well as circularization of the orbit resulting from the frictional force, is delayed. In some 
cases, the orbit stays misaligned over the disc lifetime. More massive planets and planets 
further away from the star align faster. In addition, more massive discs favour alignment. 
In section HI we discuss our results in the light of the observations that have been reported 
so far. 



2 INTERACTION BETWEEN A PLANET ON AN INCLINED ORBIT 
AND A DISC 

2.1 Gravitational potential and disc model 

We consider a planet of mass M p orbiting around a star of mass M* which is itself surrounded 
by a disc of mass Mj. The disc's midplane is in the equatorial plane of the star whereas the 
orbit of the planet is inclined with respect to this plane. We denote by / the angle between 
the orbital plane and the disc's plane. We suppose that the angular momentum of the disc 
is large compared to that of the planet's orbit so that the effect of the planet on the disc 
is negligible: the disc does not precess and its orientation is invariable (see appendix IA~|) . 
We denote by (x, y, z) the Cartesian coordinate system centred on the star and (r, <p, z) the 
associated cylindrical coordinates. The (axisymmetric) disc is in the (x, y)— plane, its inner 
radius is Ri, its outer radius is R Q and its thickness (defining the region within which the 
mass is confined) at radius r is 2H(r). 

The gravitational potential exerted by the disk at the location (r, z) of the planet is: 

.<,, ,) = —g r r r /2+ ff^; , y , (i) 

J Ri J-hJo ^/r z + r z — 2rr cos(p + {z — z) z 
where G is the gravitational constant and p is the mass density in the disc. 

We assume that p falls off exponentially with z 2 near the midplane with a cut off at 
\z\ = H(r), while decreasing as a power of r. Thus we adopt: 

p(r,z) = [eie-'VW) - l] ^L^" —fL-, (2) 

where po is a constant. We have p = at the surface of the disc, i.e. when \z\ = H, and 
p(r, 0) = po(r / R )~ n in the midplane. This ^-dependence of p has been chosen for simplicity, 
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6 J. Teyssandier, C. Terquem and J. Papaloizou 

but we note that the details of how p varies with z do not matter. What is important is the 
local disc surface density at the location the planet passes through. This largely determines 
the change in orbital energy resulting from the dynamical drag force (as is apparent from 
equation (TTTT) of section |2"U1 below) . 

The mass density given above is discontinuous at r = i?, and r = R a , as p is zero outside 
the disc. We have found that this could introduce some numerical artefacts in the calculation 
of the disc's gravitational force, so that in some runs we have replaced p by p x f(r) with: 



f(r) 



r 



10 



R„ 



20 



(3) 



The exponents 10 and 20 ensure that the edges are rather sharp, so that we get quantitatively 
the same orbital evolution whether the factor / is used or not. 

In the calculations presented below, we chose a value of the disc mass Md and calculate 
Po using Md = f J J disc pdV. For the disc's semithickness, we chose a constant aspect ratio: 



H{r) = H r, 



(4) 



where Hq is a constant. The gravitational force per unit mass exerted by the disc onto the 
planet is — V$. In appendix [Bj we give convenient expressions for — V$ in terms of elliptic 
integrals which can be readily computed. 



2.2 The Kozai mechanism 

The Kozai effect (Kozai 1962, Lidov 1962) arises when an inner body on an inclined orbit is 
perturbed by a distant companion. First derived by Kozai to study the motion of inclined 
asteroids around the Sun under perturbations from Jupiter, it has since found many appli- 
cations in astrophysics. We are interested here in cases where the inner body is a planet of 
mass M p . We denote a its semimajor axis, M' the mass of the outer companion, assumed to 
be on a circular orbit, and a' its semimajor axis. We consider the case a' 3> a. The secular 
perturbation from the outer companion causes the eccentricity e of the inner planet and the 
mutual inclination I of the two orbits to oscillate in time in antiphase provided that the 

initial inclination angle J is larger than a critical angle I c given by: 

3 

cos 2 I c = -. (5) 
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The maximum value reached by the eccentricity is then given by: 

5 \ 1/2 

1- -cos 2 Jo) , (6) 

and the time £ evo i it takes to reach e max starting from eo is (Innanen et al. 1997): 

t , / 2V 1/2 (e 



with the time r defined by: 

T= («) 3 (MAL, (8) 

\a J \M'J 2tt' V ; 

where T is the orbital period of the inner planet. If the eccentricity oscillates between e m i n 

and e max , then the period of the oscillations P osc is given by P osc = 2t evo i with e = e min in 

equation (J7J). 

Since the two orbits are well separated, the component of the angular momentum of the inner 



orbit perpendicular to the orbital plane is constant. As it is proportional to \/l — e 2 cos /, 
the oscillations of / and e are in antiphase. Hereafter, we will refer to this mechanism as the 
classical Kozai effect. 

Terquem & Ajmia (2010) found that the Kozai effect extends to the case where the inner 
orbit is perturbed by the gravitational potential of a disc, even when the orbit of the planet 
crosses the disc, provided most of the disc mass is beyond the planet's orbit. In that case, 
/ is the angle between the orbital plane and the disc's plane. They also showed that, in 
agreement with equations (I7j) and (IE]), the period of the oscillations decreased with a. It was 
also found to decrease with increasing disc mass. When the semimajor axis of the planet 
is small compared to the disc's inner radius, the evolution timescale is of the same form as 
that given by equation ((?]) but with 
R \ 3 {MA T 



where the coefficient of proportionality depends on the functional form of p and on the ratio 
Ri/ R Q . 



2.3 Friction 

When the planet crosses the disc, as it has a relative velocity with respect to the particles in 
the disc, it suffers a frictional force. There are two types of drag acting on the planet: (i) an 
aerodynamic drag, due to the fact that the planet has a finite size and suffers direct collisions 
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8 J. Teyssandier, C. Terquem and J. Papaloizou 

with the particles in the disc, and (ii) a dynamical drag, due to the fact that particles in the 
disc are gravitationally scattered by the planet. 

Adopting a drag coefficient of unity, the aerodynamic drag force per unit mass exerted 
on the planet located at (r, <p, z), or equivalently (x, y, z), can be written as: 



Taero = ~ ^Jj-H R 2 p p(r, z)v Te lV Teh (10) 

where R p is the planet's radius and v re i is the relative velocity of the planet with respect to 
the particles in the disc at the location (x,y,z). If we denote v = (v x ,v y ,v z ) the velocity of 
the planet, then v re i = (v x + yQ, v y — v z ), where Vl = a/ GM±/r 3 . 

The dynamical drag is the gravitational force exerted by the particles in the disc on the 
planet that is associated with the scattering and change of location that occurs as a result 
of the passage of the planet. The main contribution to this force comes from the particles 
located in the vicinity of the planet at the time it crosses the disc. Note that the gravitational 
force from these particles has already been included in — V$ (see section 12. ip . but as we 
have ignored the motion of the particles in the disc relative to the planet when calculating 
this force it is conservative and does not capture the change of energy of the planet's orbital 
motion. The problem is similar to dynamical friction in a collisionless medium (e.g., Binney 
& Tremaine 1987). 

When the inclination angle of the planet's orbit with respect to the disc is not very 
small, t> re i ~ a/GM^/o, with a being the semimajor axis of the planet's orbit, is supersonic. 
In this case, the dynamical friction force per unit mass acting on the planet can be written 
as (Ruderman & Spiegel 1971, Rephaeli & Salpeter 1980, Ostriker 1999): 



dyn 



-4vrG 2 M p p(r, z) 



In- 



H(r) 



Rn 



V re l 



V 



rel 



The ratio of the two frictional forces is thus: 

2 / n a- \ 2 



1 



dyn 



hi 



H(r) 



Rp 



Rp 
a 



(11) 



(12) 



where we have replaced u 2 el by GM+/a. In this paper, we will focus on planets with masses at 
least that of Neptune, for which the previous expression gives r aero /T dyn 1 for the orbital 
parameters we consider. Therefore, friction is dominated by the dynamical term even though 
a large value of unity was adopted for the drag coefficient. 
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Orbital evolution of a planet on an inclined orbit 9 

2.4 Evolution timescale 

The velocity v of the planet is, in first approximation, the Keplerian velocity around the 
star, i.e. v ~ a/ GM*/a. To simplify matters, we approximate the relative velocity by its 
vertical component so that t> re i ~ v z ~ v sin J ~ a/ GM±/a sin I. 

When the planet's orbit has a significant eccentricity, the relative velocity is in general 
larger as the planet crosses the disc closer to pericenter than apocenter and its horizontal 
components (in the disc's plane) are then important. We define a damping timescale for the 
planet in the absence of forces other than friction as: 

/ 1 dv \ 1 v . . 

T ^ = \vTt) =rv (13) 

This is the characteristic timescale on which the frictional force Tdyn damps the velocity 
of the planet. Note that Td am p — 2a(da/dt) . With the above approximation for v ie \ and 
\ln(H / R p )\ ~ 6, which corresponds to giant planets at around 10 au in a disc with H(r) = 

2.5 x lCT 2 r, as we consider below in the numerical calculations, we get: 

M 2 sin 2 / T 

Tdamp ~ 2^M p o?p{r,z)2^' ( } 

where T is the orbital period of the planet. 

To proceed further, for the purpose of getting a convenient analytical estimate of Td amp , 
we approximate the vertical dependence of the mass density by a 5-function in z and set 
p(r,z) = 25{z)H(r)p (r/R o )- n . Then we get E(r) = J H H p(r,z)dz = 2p (r/R o )- n H(r). 
Using equation (j3j), we can then calculate the disc mass to be = Yj(r)2irrdr ~ 
8ttpqHqRI/3, where we have used n = 3/2, as in the numerical simulations below. That allows 
us to express p i n term of M d , such that p ~ 3M d / (8iiH Rl). In accordance with equation 
(J2]), we identify the midplane mass density as p(r, 0) = po(r/R )~ n — 3M d (r / R )~ n / (8nH Rl) 
We now suppose that the planet's orbit is not too eccentric so that p can be evaluated at 
r = a in the expression of r damp . Finally, equation ( 114]) gives: 



' 1 ] H sm 2 I—. (15) 



aamp 9M p M d \ a J " u 2tt 
This expression may not give the correct quantitative value of the damping timescale, be- 
cause of the approximations that have been used in deriving it, but it gives the scaling of 
r damp with the different parameters. Also, although the dependence on the eccentricity has 
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not been taken into account in this expression, as pointed out above, we expect the damping 
timescale to be larger when the eccentricity is larger. 



3 NUMERICAL SIMULATIONS 

To study the evolution of the system (star, planet, disc), we use the iV-body code described 

in Papaloizou & Terquem (2001) in which we have added the gravitational and frictional 

forces exerted by the disc on the planet. 

The equation of motion for the planet is: 

d 2 v GM±r _ GM„r , , 

& = —fir ~ v$ + r — + r ^ + !V - -^jf- • (16) 

The last term on the right-hand side is the acceleration of the coordinate system based on 
the central star. It arises because the center of mass of the system does not coincide with 
that of the star. This term takes into account only the force exterted by the planet onto 
the star, and neglects the net force of the disc on the central star, which would of course 
require a calculation of the disc response to the planet's perturbation. Tides raised by the 
star in the planet and relativistic effects are included through Tt, r , but they are unimportant 
here as the planet does not approach the star closely. Equation (Tl6|) is integrated using the 
Bulirsch-Stoer method. The integrals over <ft involved in V$ are calculated using elliptic 
integrals (see appendix [B]). The integrals over r and z are calculated with the Romberg 
method (Press et al. 1993). 

The planet is set on a circular orbit at the distance r p from the star. When the planet 
does not pass through the disc, e.g. when there is no friction, the orbital energy is conserved 
and r p is equal to the planet's semimajor axis a throughout the evolution of the system. 
The initial inclination angle of the orbit with respect to the disk is Jo- In the simulations 
reported here, we have taken = 1 M , a radial power law with exponent n = 3/2 for the 
disc mass density p (see eq. [2]), a disc aspect ratio H = 2.5 x 10~ 2 and a disc outer radius 
R = 100 au. 

We will consider a planet with either (i) M p = 10~ 3 M = 1 Mj and R p = 7 x 10 9 cm 
(Jupiter), (ii) M p = 5 x 10" 5 M = 1 M N and R p = 2.5 x 10 9 cm (Neptune), (iii) M p = 
10~ 2 M = 10 Mj and R p = 8.4 x 10 9 cm, i.e. 1.2 times Jupiter's radius. These latter values 
approximately correspond to the mass and radius of the planet XO-3 b, which has been 
detected both in transit and using radial velocity measurements (Johns-Krull et al. 2008, 
Winn et al. 2008). 
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Figure 1. Eccentricity e (dotted line) and inclination angle / (in degrees, solid line) versus time (in yr) in the absence of 
friction for M p = 1 Mj, r p = 5 au, M d = 10~ 2 M , Ri = 10 au and I = 50°. 

3.1 Gravitation only 

To get a benchmark, we first consider the evolution of the planet's orbit in the case where 
friction is ignored (i.e. r aero and Tdyn are set to zero in eq. [IE]). This applies when the 
planet's distance to the star is at all times smaller than the disc inner radius (e.g., planet 
orbiting in a cavity). These calculations are similar to those performed by Terquem & Ajmia 
(2010) and a comparison is made below. 

Figure [1] shows the time evolution of the eccentricity e and inclination I for M p = 1 Mj, 
r p = 5 au, M d = 10" 2 M , R 4 = 10 au and J = 50°. 

The inclination angle varies between 7 m i n = I c and I maX ) where I c is the critical value 
of Iq below which eccentricity growth is not observed. The fact that J min = I c is illustrated 
in figure El which shows the time evolution of the inclination / for M p = 1 Mj, r p = 7 au, 
M d = 10" 2 M , Ri = 1 au and Iq varying between 23° and 50°. For these parameters, 
oscillations of the eccentricity and inclination disappear once Jo becomes smaller than 23°, 
which means that I c = 23°. For values of Iq well above I c , the inclination is seen to oscillate 
between I min = I c and / max - When I decreases however, the amplitude of the oscillations is 
such that I min gets a bit larger than I c . 
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Figure 2. Inclination / (in degrees) versus time (in yr) in the absence of friction for M p = 1 Mj, r p = 7 au, = 10~ 2 Mq, 
Ri = 1 au and Io = 50° (solid line), 40° (dotted line), 30° (short-dashed line) and 23° (long-dashed line). For these parameters, 
Ic = 23°. No oscillations are present for Iq < I c . 

As already noted by Terquem & Ajmia (2010), since I c depends on the gravitational force 
exerted onto the planet, it varies with the initial position of the planet r p . When r p <C Ri, 
the conditions are similar to the classical Kozai effect and the gravitational force from the 
disc can be approximated by a quadrupole term. In this case, it can be shown that I c = 39° 
(eq. [5]). When r p is larger though, the quadrupole approximation is not valid anymore and 
I c gets smaller. This is illustrated in figure [3J which shows the critical angle I c as a function 
of r p , which here is the same as the planet's semimajor axis since the orbit is initially 
circular and there is no energy dissipation. These calculations are done for M p = 1 Mj, 
M d = 10~ 2 M , Ri = 10 au and r p varying between 1 and 35 au. For these parameters, the 
Kozai effect disappears when r p becomes larger than ~40 au, as most of the mass is then 
no longer outside the orbit of the planet (Terquem & Ajmia 2010). 

For comparison, we have also run the case displayed in figure [1] by modelling the disc 
in two dimensions only, as in Terquem & Ajmia (2010). The gravitational potential is then 
calculated as arising from a distribution of mass with surface density (0 ^ r ~ v ■ We take 
p = 1/2 as in the 3D calculations we have S(r) ~ p{r)H{r) oc r -1 / 2 . There is very good 
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Figure 3. I c (in degrees) versus r p (in au) in the absence of friction for M p = 1 Mj, = 10 — 2 Mq, Ri = 10 au and r p 
varying between 1 and 35 au. For r p <C Ri au we recover the classical Kozai value I c = 39° . 

agreement between the two sets of calculations. The extrema of I and e and the period of 
the oscillations are roughly the same. 

In the 3D case, we have also checked that the vertical structure of the disc does not affect 
the oscillations. Indeed, when the disc mass is kept constant but H is varied in equation @, 
the oscillations are unchanged. 

3.2 The effect of friction 

We are now going to study the effect of friction on the evolution of the planet's orbit by 
taking into account r aero and Tdyn in equation ( floT) . 

Figure H] shows the time evolution of the orbital parameters /, e, semimajor axis a and 
distance to pericenter a(l — e) for M p = 1 M N , r p = 7 au, M d = 1CT 2 M , = 1 au 
and for two different values of I Q . For these parameters, I c = 23° (the critical angle does 
not depend on the mass of the planet an therefore is the same as in fig. [2]). As discussed in 
section I2.4[ the damping timescale increases with / and e. Therefore, when I a > I c , as at 
all times one of these parameters has a large value because of the Kozai cycle, the damping 
timescale is longer than in the case I Q < I c and e = 0. This is illustrated in figure HI where 
we see that a decreases much more rapidly when I = 20° than when I Q = 50°. In the 

© 2012 RAS, MNRAS OOO.ITH281 



14 J. Teyssandier, C. Terquem and J. Papaloizou 

case where there is no Kozai cycle, the orbit stays circular and I is damped faster as it 
becomes smaller, so that the orbit aligns with the disc on a very short timescale. When 
I Q > I c , Kozai oscillations are present as expected, but the oscillations are damped because 
of friction. We observe that damping is much less efficient than in the I Q = 20° case, even 
though the inclination / does reach values that are not much larger than 20° and for which 
the damping timescale would be similar if the orbit were circular. When the Kozai cycle is 
present though, the eccentricity of the orbit is large when I is minimum, which results in 
the damping timescale being maintained at large values. Contrary to what happens when 
I Q < I c then, the damping timescale stays roughly constant. This is illustrated in the plot 
of figure H] that shows a versus time. Ultimately, if the disc were present long enough, the 
orbit of the planet would align with the disc and would be circularized (I and e vanish). For 
a Neptune mass planet and the parameters used here though, we find that misalignment 
can be maintained over the the disc lifetime, which is of a few million years. A Jupiter mass 
planet would align faster, as discussed below. 

These calculations show that, by pumping up the eccentricity of the planet's orbit and 
maintaining either / or e at large values, the Kozai effect provides a way of delaying alignment 
with the disc and circularization of the orbit. 

We see in figure S] that the period of the oscillations decreases with time. If there were 
no dissipation, this period would stay constant. In the classical Kozai cycle, when the orbit 
of the inner planet has a semimajor axis a, the time it takes to reach e max starting from 
e min , is proportional to a~ 3//2 In (e max /e min ) (see eq. [7] and [8]). Although this formula has 
been derived for a quadrupolar gravitational potential, it may be expected to give a general 
trend in the disc case as well. Friction reduces a, which tends to increase this timescale. But 
the amplitude of the cycle also decreases (e min increases and e max decreases), so that the net 
effect is a shorter timescale. 

As pointed out in section I2"3j the expression for dynamical friction given by equation ( TTTj) 
is valid only when the relative velocity of the planet with respect to the particles in the disc 
is supersonic. This is not the case when I is close to zero and, therefore, the part of the 
curves in figure H] corresponding to I close to zero can be interpreted only qualitatively. 
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Figure 4. Inclination / (in degrees, upper left plot), eccentricity e (upper right plot), semimajor axis a (in au, lower left plot) 
and pericenter distance a(l — e) (in au, lower right plot) versus time (in yr) for M p = 1 Mn, r p = 7 au, = 10 — 2 Mq, 
Ri = 1 au and Iq = 50° (solid line) and 20° (dotted line). For these parameters, I c = 23°. When Iq = 20°, e = throughout 
the evolution so the curve does not show on the plot. The Kozai effect provides a way of pumping up e and maintaining either 
e or I at large values, and therefore delays circularization and alignment with the disc of the orbit. 

3.2.1 Influence of the planet's mass: 

We now study the influence of the planet's mass on the evolution of the orbit. Figure 
shows the time evolution of the inclination /, eccentricity e and semimajor axis a for the 
same parameters as in figure HI i.e. r p = 7 au, Md = 1CT 2 M , Ri = I au, I Q = 50° and three 
different values of M p . As I c = 23° for these parameters, we are in the regime where Kozai 
cycles are present. As shown by equation ( 1T5"|) . the damping timescale r damp is proportional 
to 1/M P . The period P osc of the Kozai cycle, on the other hand, does not depend on the 
planet's mass (see eq. [S] and [5]). We therefore expect friction to be very efficient for high 
mass planets for which the damping timescale rd am p would be smaller than P sc, and much 
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less efficient for low mass planets for which Td am p 3> -Pose- This is borne out by the results 
displayed in figure EJ Friction dominates the evolution for the planet with M p = 10 Mj, 
for which 7d amp is smaller than P osc . For a Jupiter mass planet, the two timescales are 
comparable and the orbit aligns after a few oscillations, i.e. after a few 10 5 years. For a 
Neptune mass planet, rd amp P osc and / decreases very slowly. In that case, the orbit 
would stay misaligned over the disc lifetime. 

As stated above, P osc does not depend on the planet's mass. However, we see from figure E] 
that the oscillations for M p = 1 Mj and M p = 1 Mn do not have the same period. This is 
because, as noted above, P osc decreases as damping reduces the amplitude of the oscillations. 

The fact that dynamical drag leads to faster alignment of the orbits for more massive 
planets was seen in the numerical simulations of Rein (2012) who considered planets on 
highly inclined orbits but without taking into account the gravitational interaction with the 
disc. 

3.2.2 Influence of the disc's mass: 

Figure [6] shows the time evolution of the inclination / and semimajor axis a for M p = 1 Mn, 
r p = 7 au, Ri = 1 au, I = 50° (as in fig. @| and three different values of M d ranging from 
10 -2 to 0.1 M . Terquem & Ajmia (2010) showed that I c does not depend on when 
the other parameters are kept fixed, so that we have I c = 23° here. Since Jo = 50°, we are 
therefore in the regime of Kozai cycles. From equation ( fl5l) . we see that rd amp oc M^j~ . On 
the other hand, we also have P osc cx , as expected from the classical Kozai effect and 
verified in the disc's case by Terquem & Ajmia (2010). It follows that alignment of the orbit 
with the disc, when it occurs, should happen after the same number of oscillations whatever 
the disc's mass. This is confirmed by the curves in figure [6] which are all a time-scaled 
version of each other, with the more massive disc showing the shortest evolution timescale. 
It is therefore easier to keep a planet on an inclined orbit in the presence of a less massive 
disc. 

3.2.3 Influence of the planet's initial semimajor axis: 

Finally, we study the effect of the initial semimajor axis on the orbital evolution. Figure [7] 
shows the time evolution of the inclination I for M p = 1 Mj, Md = 10~ 2 M , Ri = 10 au, 
Io = 60° and four different values of r p ranging from 5 to 20 au. For these parameters, we 
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Figure 5. Inclination / (in degrees, upper plot), eccentricity (middle plot) and semimajor axis a (in au, lower plot) versus 
time (in yr) for r p = 7 au, = 10 — 2 Mg, Ri = 1 au, Io = 50° (same parameters as in fig- and M p = 1 Mn (solid line), 
1 Mj (dotted line) and 10 Mj (dashed line). For these parameters, I c = 23° < Iq so that we are in the regime where Kozai 
cycles are present. Friction damps the oscillations more efficiently for larger mass planets. Alignment of the orbit with the disc 
may not happen over the disc lifetime for smaller mass planets. 



are in the regime of Kozai cycles. For r p = 5 au, the planet stays in the disc's inner cavity so 
that it never experiences frictional forces when it goes through the disc's plane. The orbital 
evolution is therefore that of a Kozai cycle with no friction and a minimum angle of 37°, 
consistent with the results of figure [31 For r p = 8 au, although the distance to the star does 
get larger than Ri when e > 0.25, as it happens the planet is always in the disc's inner cavity 
when it crosses the disc's plane, so that here again friction does not play a role. Like in the 
previous case, the minimum angle for this cycle is consistent with the results of figure [31 As 
noted above, P osc oc a -3 / 2 in the classical Kozai case. We verify here that P osc does indeed 
decrease when a increases (in agreement with Terquem & Ajmia 2010). For r p = 12 and 
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Figure 6. Inclination I (in degrees, upper plot) and semimajor axis a (in au, lower plot) versus time (in yr) for M v = 1 Mn, 
r p = 7 au, Ri = 1 au, To = 50° (same parameters as in fig. [4j| and Md = 0.1 Mq (dashed line), 5 X 10~ 2 Mq (dotted line) and 
10 — 2 Mq (solid line). For these parameters, I c = 23° < Iq so that we are in the regime where Kozai cycles are present. It is 
easier to keep a planet on an inclined orbit in the presence of a less massive disc. 



20 au, the planet crosses the disc so that the oscillations are damped by the frictional force. 
As shown by equation ffT5l) . the damping timescale Td am p does not depend on a. This is 
consistent with the curves displayed in figure [7] at early times, before the eccentricity of the 
orbits grow, as they show that the damping timescale is indeed the same for r p = 12 and 
20 au. But as the inclination gets smaller for the larger value of r p (I m i n is smaller and the 
planet spends more time in the disc), 7d amp gets also shorter for this value of r p . Alignment 
of the orbit with the disc is then faster at larger distances from the star. 
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Figure 7. Inclination / (in degrees) versus time (in yr) for M p = 1 Mj, = 10 — 2 Mq, Ri = 10 au, Io = 60° and r p = 20 au 
(solid line), 12 au (dotted line), 8 au (short dashed line) and 5 au (long dashed line). For these parameters, we are in the regime 
of Kozai cycles. Alignment of the orbit with the disc is faster at larger distances from the star. 

4 DISCUSSION 

We have investigated the dynamics of a planet on an orbit inclined with respect to a disc. If 
the initial inclination of the orbit is larger than some critical value, the gravitational force 
exerted by the disc on the planet leads to a Kozai cycle in which the eccentricity of the orbit 
is pumped up to large values and oscillates with time in antiphase with the inclination. 
On the other hand, when the planet goes through the disc, it suffers a frictional force that 
results in a loss of orbital energy. As a consequence, the inclination and the eccentricity of 
the orbit are damped and the semimajor axis decreases. The goal of this paper was to study 
on what timescale orbits inclined with respect to a disc would align with it. 

The calculations presented in this paper show that, by pumping up the eccentricity of 
the planet's orbit and maintaining either I or e at large values, planets in orbits undergoing 
Kozai cycles maintain large velocities relative to the disc as they pass through it, so delaying 
alignment with the disc and circularization of the orbit. 

For the parameters used in this paper, which are typical of protostellar discs, it was found 
that Neptune mass planets would remain on inclined orbits over the disc lifetime. Jupiter 
mass planets, however, tend to align faster, as the damping timescale is shorter for more 
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massive planets. Note however that we have not taken into account the fact that the disc 
dissipates progressively over time. As damping is less efficient in less massive discs, Jupiter 
mass planets could remain misaligned if the disc's mass were decreasing sufficiently fast. We 
have also found that alignment of the orbits was faster at larger distances from the star. 

So far, the only planets that have been found on inclined orbits are rather massive 
(with a mass ~ Mj) and on short period orbits. This of course is a result of observational 
bias, as the inclination is measured so far only for transiting planets. As the results of this 
paper suggest, if these planets had been on inclined orbits when the disc was present, they 
probably would have aligned if they had crossed the disc. Therefore, either (i) the orbits 
became inclined after the disc had dissipated, (ii) the disc in which the planets formed was 
misaligned with the stellar equatorial plane, or (iii) the planets formed on inclined orbits 
with short enough periods that they crossed the disc's plane only in an inner cavity. Jupiter 
mass planets formed on inclined orbits at large distances from the star would be expected 
to have aligned with the disc unless the formation took place near the end of the life of the 
disc. 

The planets in the Kepler sample seem to have both low inclinations (Fabrycky et al. 
2012) and low eccentricities (Kane et al. 2012). Radial velocity surveys also show lower 
eccentricities for lower mass planets. According to the results of our paper, this strongly 
suggests that these planets, with masses mainly at most ~ Mpj, have formed and stayed in 
the original disc. Had they been in a sufficiently inclined orbit at some point, this would have 
recurred over the disc lifetime, and the orbit would have had episodes of high eccentricity. 
Thus a population in both highly eccentric and inclined orbits would be expected. Measure- 
ment of inclination angles for longer period orbits will enable the evaluation of proposed 
aspects of planet formation scenarios. 

Finally, we remark that in the present paper we have considered only one planet inter- 
acting with the disc. In a subsequent paper, we will investigate the dynamics of multiple 
systems. 
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APPENDIX A: RESPONSE OF THE DISC TO THE GRAVITATIONAL 
INTERACTION WITH THE PLANET 

We here estimate the warping response of a disc of the type we consider to a planet in an 
inclined circular orbit. The disc is assumed to contain significantly more angular momentum 
than the planet and obey a barotropic equation of state. We show that provided the inverse 
of the orbital precession frequency is less than the local disc sound crossing time and the 
mass of the planet is less than the disc mass in its neighbourhood, the range of inclinations 
excited is expected to be small. 
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Al Governing equations 

The basic governing equations are the equations of continuity and motion for a barotropic 
gas in the form 



^ + V-pv = 0, (Al) 
<9v 1 

+ v-Vv = — VP-V$-V$*, (A2) 
at p 



where = — GM*/\r\ is the gravitational potential due to the central star and 



* = GM P 

y/r 2 + R 2 - 2rRcos(4> - <j) p ) + (z - z p ) 2 



is the potential due to the planet which is treated as a perturbation for which we calculate 
the linear response below (the indirect term does not contribute to the warping and so may 
be dropped). Here the cylindrical coordinates of the planet are (R,<p p ,z p ). 
We adopt a Fourier decomposition in azimuth and time of the form 

$ = ®™ exp [i(m0 + u Ptm t)] + cc, ( A4) 

m>0 

where +cc indicates the addition of the complex conjugate, m is the azimuthal mode number 
and U) p ^ m is an associated frequency corresponding to a pattern speed —uj p ^ m /m. For global 
warps we are interested in m — 1 and pattern speeds that correspond to a slow precession 
of the planetary orbit. The precession period of the line of nodes, which for the purposes 
of this section is assumed to be given, is related and comparable to the period of Kozai 
oscillations when these occur. We adopt the time or orbit average of the coefficient $i 
which is appropriate for a discussion of the secular evolution of global warps. We perform 
a corresponding Fourier decomposition for the response perturbations to the disc which are 
taken to have a <p and t dependence through a multiplicative factor exp[i(m<f)+L) P)m t)]. From 
now on this factor will be taken as read and we drop the subscript m on quantities as this 
is taken to be unity. 
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Figure Al. The left hand panel shows the inclination, normalized by the factor 10 3 (M p / 'M*)(0.05R/ H (R)) 2 uj p R 3 ^ 2 /\/G M», 
as a function of r/R for the calculation described in the text. This was calculated under the condition that it tended to zero as 
r — > oo. The right hand panel shows the function T(r / R)/T(0) which represents the magnitude of the torque due to the planet 
acting on the disc in the radial interval (r, oo) normalized by its value as r — > 0. For both panels the inclination of the circular 
planetary orbit of radius R to the plane of the disc was it/ 4. 



A2 The disc inclination response 

Denoting perturbations with a prime, linearization of the equations of motion (IA2I) for the 
response to the perturbing potential, <3>, gives 

i(u p + tt)v' r - 2ttv' = 



dW 



i(u p + + — v' r 
i(u p + VL)v' z 



dr 
\W_ 

r 

dW 

dz 



(A5) 



where W = P' / p + $ and n is the epicyclic frequency. Solving for the velocity perturbations, 
we obtain 

_. (up + n)dW/dr + 2VLW/r 

Vr ~ _1 k 2 - (up + ny 

(K 2 /(2Sl))dW/dr + (u p + Q)W/r 

V * ~ - {Up + W { } 

As we are interested in a disc that is close to Keplerian rotation, with u p Q, and k ~ Q, 
we neglect u p and set k = Q in the numerators above, and replace k 2 — (u p + O) 2 by 
2Q(k — Up — Q) in the denominators to obtain 
. dW/dr + 2W/r 

Vr ~ _1 2 ( K - up - n) 

dW/dr + 2W/r 

v <$> ~ T7 rvT ' 

v 4{K — Up — U) 
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Following our previous work (eg Papaloizou & Terquem 1995, Larwood et al. 1996, Nelson &; 
Papaloizou 1999) we seek a solution for which v' z is independent of z to within a correction 
of order (H/r) 2 . Then to within the same order of accuracy we may integrate the linearized 
z component of the equation of motion to give 

W = -i(oo p + Q)zv' z . (A8) 
We now write down the linearized continuity equation in the form 

iK + ny-» ) = 

c. 

where we have used P' = c 2 p', with c 2 = dP/dp. Multiplying ( 1A9I) by z and integrating over 
the vertical extent of the disc, we get 

i(u p + n)$pz , f 00 (u p + Q) 2 v'pz 2 , f°° , , „ / f°° . ,\ 

— - ; - dz = / ^ ! 2 z — dz- j v' z pdz+ V_l-(/ -zv^pdz J , (AlO) 

where the perpendicular velocity perturbation is = (v' r ,v'^,0). Making use of equations 

(1A7I) . (IA8p with w p neglected in the latter, and vertical hydrostatic equilibrium for the 

unperturbed state, we obtain 

d ( pVL dg\ , Auj p Y>g 2ir 2 f°° 9$ J 2ir 2 S ^9$ 



arV^ + o-Kar; + n gmJ^ dz dz ~ gm* [dz ) z=Q } (A11) 

where 



//= / p-2 2 dz, (A12) 

g = r 2 Qv z /(r 3 Q 2 ) = ^/(rfl) and terms of order a; 2 have been neglected. Note that con- 
sistently with the approximations used here, where possible, we have adopted Keplerian 
rotation and taken r 3 Q 2 = GM* to be a constant in this expression. Thus 2g becomes the 
local inclination of the disc (the factor of two arises from the form of the Fourier decompo- 
sition f lA4j) ). We remark that (I A 1 1 1) may also be written as 



d ( pQ dg\ , Au^g 1 dT 



dr \tu p + n — ndr J Q nGM* dr 
where, for orbits with line of nodes coinciding with the y axis as considered below, T is the 
torque due to the planet acting on the disc in the radial interval (r, oo) in the x direction. 
Thus T — > as r — > oo. 

A3 Solution for the disc inclination 

We wish to solve ( 1A13I) for the inclination of the disc induced by a planet on an inclined 
orbit. As the unforced problem has solutions consisting of long wavelength bending waves 
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propagating with a speed that is a multiple of the sound speed (Nelson & Papaloizou 1999), a 
complete solution requires knowledge of the complete structure of the disc including bound- 
ary details which are beyond the scope of the modelling in this paper. To deal with this 
situation we assume the disc extends to large radii and has a much larger angular momen- 
tum content than the planet. We then look for a solution for which the inclination response 
is localized away from large radii having g — > 0, and dg/dr — > 0, for r — > oo. These solutions 
are in fact regular as r — > for the power law discs we adopt here. However, inner boundary 
effects could possibly cause the excitation of a freely propagating bending wave that should 
be added. But this should not affect an estimate of the scale of the warping. To simplify 
matters further, we shall take k = Q as expected for a constant aspect ratio disc, under the 
gravitational potential due to a central point mass, and assume that u p /Q is much less than 
H/r. The latter assumption, which is expected to lead to mild warping (eg Larwood et al. 
1996 and see below) enables us to neglect the term oc g in equation (IA13j) which can then 
be easily integrated to give 
yufi dg 1 

~Up~dr~ = 7rGM* ( } 

and accordingly 

9 = -^] r ^MM/ r - (A15) 

We make a rough estimate of g as determined by (1A15|) by setting /i = Eif 2 and T ~ 
irGMpJ^r ~ M p r 2 u p fl, where r is evaluated at a location where the planet intersects the 
disc. Then we estimate 

g „ „ (A16) 

This implies that the inclination range is in general small. The first factor in brackets mea- 
sures the square of the product of the precession frequency and local sound crossing time 
which is expected to be less than unity and so lead to only small warping (see Papaloizou & 
Terquem 1995, Larwood et al. 1996, Nelson & Papaloizou 1999). The second factor in brack- 
ets is expected to be the ratio of the planet mass to the characteristic disc mass contained 
within a length scale comparable to that of the orbit, which is expected to be less than of 
order unity particularly for low mass planets. 

We have evaluated the solution given by (1A15j) for a planet in a circular orbit of radius 
R, inclined at 45 degrees to the plane of the disc and with line of nodes coinciding with the 
y axis. We took S oc r" 1 / 2 as in the main text and as the solution scales with the inverse 
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square of the disc aspect ratio, H/r, that is left as a parameter. The left hand panel of 
Fig. ED shows the inclination in units of 10 3 (M p / 'M*)(0.05i2/ 'H{R)) 2 u p R^ 2 / 'y/UM*, as a 
function of r/R. The maximum value of this is of order unity indicating that, even for a 
Jovian mass planet in a disc with aspect ratio 0.05, the characteristic value of the inclination 
range is of order co p R 3 ^ 2 / y/GM* which is expected to be of the order of the ratio of the mass 
of the disc to that of the central star and thus a small quantity. The right hand panel of Fig. 
IA1I shows the cumulative torque T(r/R)/T(Q) measured in the sense of increasing inwards. 
This indicates that the torque falls off rapidly at large radii. 



APPENDIX B: EXPRESSION OF THE GRAVITATIONAL FORCE DUE 
TO THE DISC IN TERMS OF ELLIPTIC INTEGRALS 

Here we develop expressions for the gravitational force per unit mass exerted by a disc on a 
planet that may be passing through it in terms of elliptic integrals. The resulting expressions 
require the evaluation of two dimensional integrals with integrands that at worst contain a 
logarithmic singularity which is readily managable numerically. 

The gravitational potential exerted by the disk at the location (r, z) of the planet is given 
by equation (fTJ) as 



J R t J-hJo yr + r z — Zrr cos <p' + [z — zy 
As the disc is axisymmetric, the gravitational force per unit mass exerted by the disc on 
the planet has only a radial and a vertical components, given by —d^/dr and —d^/dz, 
respectively. To calculate these, we first note that $ is solution of Poisson's equation, which 
is given by 

d 2 $ 13$ d 2 $ m 2 $ A „ . s 
or z r or oz r z 

with m = 0. In the general nonaxisymmetric case, when the ip dependence of p is through 
a factor exp(i<£>), the azimuthal number m is nonzero and we have d 2 Q/dip 2 = — m 2 $. We 
now differentiate equation ( 1B2|) . with m = 0, with respect to r to obtain to 



d 2 + ld_ fm>\ _ + d*_ _ AnG dp , B3 , 

dr 2 \dr J r dr \ dr J r 2 dr dz 2 \dr J dr 

This shows that <9$/<9r satisfies Poisson's equation with m = 1 and dp/dr as source term. 
The solution can be written as 
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d $ ^ i>R [H ^ &( r ',z')r'dr' cos cp'dz'dcp' 



— = -G / / / - drV = = + Bi + B . (B4) 

dr Jr. J_ h J v / r 2 + r /2 _ 2rr > cos 0' + (z - z') 2 

Here the domain of integration is the interior of the domain containing the mass distribution 
and Bi and B a are boundary terms that have to be taken into account when the disc's mass 
density is given by equation (TSJ), as in that case dp/dr is infinite at the radial boundaries of 
the disc. When the disc's mass density is made continuous by multiplication by the factor 
/ defined in equation ([3]), these boundary terms are not included. 
Expression (1B4j) can be recast in the form: 



<9$ 
dr 
with 



^ = -G [ ° [ |V, z')H r (r, r', z - z')r'dr'dz' + B t + B , (B5) 
or Jr. J^ h dr 



- f -m * 7™« + < m (B6) 

Jo v r + r — ^ rr cos <P + {z — z ) 
We define a 2 = r 2 + r' 2 + (z — z') 2 , b 2 = 2rr' and u = 2b 2 / (a 2 + b 2 ). It is then straightforward 

to show that: 



H r (r, r', z — z') 



4Va 2 + b 2 \ a 2 

;K(u) — E(u) 



b 2 



a 2 + b 2 ' 



(B7) 



where K and E are the elliptic integrals of the first and second kind, defined as: 



K{m) = / - (B8) 
Jo v 1 — m sin 9 



tt/2 



E(m) = / Vl-msin 2 ^, (B9) 
Jo 

with m < 1. 

We calculate the boundary terms B { and B D by assuming that p is continuous and supposing 
that p increases from to p(Ri, z) over a distance Ar — » at the inner edge, and decreases 
from p(R Q , z) to over the same distance at the outer edge, i.e.: 

B,i = -G ( [ ^-(r',z')H r (r,r',z- z')r'dr'dz', (BIO) 



Ri-Ar J-H or 

-G I p(R h z')RiH r (r,Ri,z- z')dz . (Bll) 
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Similarly: 



B 



~G / / ^(r',z')H r (r,r',z-z')r'dr'dz', (B12) 

JRo J-h or 

l>H(R ) 

= +G p(R ,z')R H r (r,R ,z- z')dz' . (B13) 

J-H(R ) 

We calculate dQ/dz in a similar way by differentiating Poisson's equation (1B2[) . with 
m = 0, with respect to z, which leads to: 



<9r 2 \dz J r dr \ dz J dz 2 \dz J dz 

This shows that d^/dz satisfies Poisson's equation with m = and dp/dz as source term. 
The solution can be written as 



<9$ [Ro [H ,2* ^(r',z')r'dr'd(b'dz' 

— = -G / / / fel ! ; y . B15 

<9z 7^ 7_ hJo y/r 2 + r /2 - 2rr' cos ^ + (i - z') 2 

This expression can be recast in the form 



d$ [ Ro [ H dp 

-j— = — G / / —(r ,z)H z (r,r ,z — z)r dr dz , (Bf6) 

oz JRi J-H oz 

with: 



H z (r,r',z-z')= f** S= (B17) 

Jo \/ r 2 + r /2 — 2rr' cos 0' + (z — z') 2 

or, in term of the elliptic function K: 



H z (r,r',z-z') = (B18) 
Va 2 + o 2 
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